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Abstract 



Lattice Boltmzmann Methods (LBM) have been proved to be very effective methods for computational aeroacous- 
tics (CAA), which have been used to capture the dynamics of weak acoustic fluctuations. In this paper, we propose 
a strategy to reduce the dispersive and disspative errors of the two-dimensional (2D) multi-relaxation-time lattice 
Boltzmann method (MRT-LBM). By presenting an efl'ective algorithm, we obtain a uniform form of the linearized 
^SJ Navier-Stokes equations corresponding to the MRT-LBM in wave-number space. Using the matrix perturbation the- 

ory and the equivalent modified equation approach for finite diff'erence methods, we propose a class of minimization 
problems to optimize the free-parameters in the MRT-LBM. We obtain this way a dispersion-relation-preserving 
O LBM (DRP-LBM) to circumvent the minimized dispersion error of the MRT-LBM. The dissipation relation preci- 

sion is also improved. And the stability of the MRT-LBM with the small bulk viscosity is guaranteed. Von Neuman 
^ analysis of the linearized MRT-LBM is performed to validate the optimized dispersion/dissipation relations consid- 

S ering monochromatic wave solutions. Meanwhile, dispersion and dissipation errors of the optimized MRT-LBM are 

O quantitatively compared with the original MRT-LBM . Finally, some numerical simulations are carried out to assess 

. the new optimized MRT-LBM schemes. 

. Sh Keywords: Computational aeroacoustics. Lattice Boltzmann, DRP-LBM, Dispersion, Dissipation, Von Neumann 

C/^ analysis 

>. 

1. Introduction 

J^ The lattice Boltzmann method (LBM) has emerged as a very eff'ective methodology for the computational model- 

_i- ing of a wide variety of complex fluid flows [1]. Recent researches dealing with dispersion and dissipation relations 

\f^ of the LBM have shown that the LBM possesses the required accuracy to capture weak acoustic pressure fluctuations 

"^ [2, 3, 4, 5]. The analysis indicated that the simple LBS possess lower numerical dissipation than the aeroacoustic op- 

h"^ timized schemes of high-order schemes for Navier-Stokes equations [3]. The second-order accurate LBM has better 

^^ dispersion capabilities than the classical Navier-Stokes schemes with 2nd-order accuracy in space and 3 -step Runge- 

^-H Kutta in time [3]. However, the dispersion error in the LBM is higher than that in the finite diff'erence method with 

^^ 3rd order spatial discretization and the 4th order time discretization, and also higher than that in dispersion relation 

^ preserving (DRP) 6th order accurate schemes. From the view of the numerical computations, for a given dispersion 

error, the LBM is faster than the high-order schemes for Navier-Stokes equations [3]. It has been shown that the dis- 
persion error can be considered as a weakness of the LBM. For the classical lattice Boltzmann model (BGK-LBM), 
C^ it is impossible to reduce the dispersion/dissipation errors. Meanwhile, it is reported that the original MRT-LBM 

[6, 7] and the BGK-LBM have exactly the same dispersion error and there exists a high dissipation of the acoustic 
modes for the MRT-LBM [3]. Here, the original MRT-LBM means we use the relaxation parameters recommended 
by Lallemand and Luo [6] . Furthermore, because of a high value of the bulk viscosity, the original MRT-LBM has a 
better stability compared with the BGK-LBM [3]. So, if the dissipation error of the MRT-LBM can be reduced and 
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the stability can be guaranteed, it will be a good choice to simulate the acoustic problems using the MRT-LBM. The 
MRT-LBM can be optimized thanks to the existence of free parameters. By means of a Taylor series expansion [8, 9], 
it is easy to establish the relation between the linearized MRT-LBM (L-MRT-LBM) and the linearized Navier-Stokes 
equations (L-NSE) with the high-order truncation. The relations between the L-MRT-LBM and the L-NSE offer us 
a way to detect the influence of free parameters on the dispersion/dissipation relations. In the limit of linear acoustic, 
von Neumann analysis is a reliable tool to recover the dispersion and dissipation relations of the LBS. It is noted that 
this famous analysis method has been revisited and extended [3, 10]. Considering plane wave solutions, the relation 
between the wave-number k and the wave pulsation co is described numerically. Then, the influence of free parameters 
on acoustic modes and shear modes will be analyzed. We propose a class of optimization strategies to minimize the 
dispersion/dissipation errors based on the matrix perturbation theory and the modified equation approach, leading to 
the definition of dispersion/dissipation-relation-preserving MRT-LBM (D2RP-LBM) schemes. 

The basic idea of the DRP-LBM is significantly diff'erent from the idea of the classical DRP-schemes correspond- 
ing to finite difference schemes [11]. The classical DRP-schemes were established considering the finite difference 
approximation of the first derivative dfjdx at the node of an uniform grid in wave-number space [11]. The classical 
DRP finite difference approximation of the first derivative df/dx in wave-number space is given by 
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The effective wave-number of Eq. (1) can be rewritten as follows 

_. M 
j=-N 



In the physical spaces, the expression of Eq. (1) is given as follows 






In order to minimize the dispersion error, the following integral error E is defined [11] 
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\ab.x - QfAxp (i(aAx). (4) 



The classical DRP schemes only focus on proposing a best approximation of the first-derivative df jdx on an uniform 
mesh. In this paper, the proposed method to minimize the dispersion/dissipation error focus on obtaining a best global 
approximation of the exact L-NSE systems based on the recovered L-NSE by the L-MRT-LBM in wave-number 
spaces. Consequently, the resulting approximation addresses the best approaching relation between the exact L-NSE 
and the recovered L-NSE, but does not address individual derivatives. Formally, the research in this paper is dedicated 
to establishing the approximation between the following equation systems in wave-number spaces 

dtW = B'W (exact L-NSE), dtW = 5,„ • W + 0{6f) (recovered L-NSE), (5) 

where the square matrices B and En are functions of wave-number vector k, where Bn can be regarded as a per- 
turbation of B. The vector W denotes the perturbed macroscopic fluid flow quantities (density, momentum) in the 
linearized-NSE systems. In the paper, in order to reduce the dispersion error of the MRT-LBM for zero-mean flow, 
the truncation error is up to 0{dt^). Meanwhile, in order to reduce the dissipation error of the MRT-LBM for uniform 
flows, the truncation error is up to Oidt"^). 

Furthermore, the proposed derivation of higher-order Taylor expansions of the MRT-LBM is very lengthy and 
complicated [8]. However, this derivation still pave a new way for establishing the relation between the L-MRT-LBM 
and the L-NSE. In order to avoid using this complex derivation, a new more effective and easy-to-use recursive algo- 
rithm is proposed to recover the L-NSE by the L-MRT-LBM. This recursive algorithm is established in wavenumber 



space. By this algorithm, the corresponding Hnearized macroscopic equations are expressed by an easy-to-handle 
matrix form. The optimization strategy is precisely built on the basis of this matrix equation. Finally, von Neumann 
analysis and numerical tests are implemented to assess the optimized MRT-LBM. 

In next section, the methodology used to establish the transformation relation from the L-MRT-LBM to the L- 
NSE is given. The optimization strategies corresponding to the matrix equation are studied in Section 3. In Section 4, 
by von Neumann analysis, the DRP-LBM schemes are analyzed in spectral spaces. In the last section, the optimized 
MRT-LBM schemes are validated by benchmark problems. 

2. Methodology for bridging from the MRT-LBM to the hnearized Navier-Stokes equations 

In this section, the basic theory of lattice Boltzmann schemes briefly reminded. Then, the linearized LBM is 
introduced and a method to establish the relation between the linearized LBS and the L-NSE is proposed. 

2.1. Lattice Boltzmann schemes 

The evolution equations of distribution functions of the basic lattice Boltzmann schemes are written as follows 

fi{x + Vidt, t + 6t) = fi(x, t) + Kij [fp\x, t) - fj(x, O) , < /, j < N, (6) 

where v/ belongs to the discrete velocity set ^, fi(x, t) is the discrete single particle distribution function corresponding 
to Vi and ff^ denotes the discrete single particle equilibrium distribution function. 6t denotes the time step and N +\ 
is the number of discrete velocities. A/y is the relation matrix. From here, the repeated index indicates the Einstein 
summation for to N except for the special indications. Let X e M^ (J denotes the spatial dimension) denote the 
lattice system, the following condition is required [8] 

X + Vjdt G £, (7) 

that is to say, if x is a node of the lattice, x + Vj6t is necessarily another node of the lattice. 
For the BGK-LBM, the relaxation matrix is set as follows 

A.7 = ^^.7, (8) 

where s is related to the relaxation frequency of the BGK-LBM. 
The standard MRT-LBM has the following form [8, 6] 

nii = Wi = mf'^\0<i<d, (9) 

and 



m. 



(x -h StVj, t + St) = mi(x, t) -h Si {mf'^\x, t) - mi(x, t)),d+l <i<N, (10) 

where the index / doesn't indicate the summation. According to the work of Lallemand and Luo [6], the relaxation 
parameters in Eq. (10) must satisfy the following stability constrain 

^,G (0,2), J -hi </<A^. (11) 

The relaxation matrix A associated with Eqs. (9) and (10) is defined by 

A = M-^^M, (12) 

where 5" is a diagonal matrix which is related to the relaxation parameters of the MRT-LBM. M = (M,/) is 

the transformation matrix, which satisfies the following basic conditions [6] 

Moj = 1, M^j = vj, (1 < a < d). (13) 

The macroscopic quantities are defined by [6, 8] 

mt = Mtjfj, mf''^ = Mijfp\ (14) 



22. The Linearized MRT-LBM and the higher-order linearized NSE 

The equilibrium function mf^ (J + 1 < i < N) is the function of conservative quantities Wi (Wt = rnf^\ < i < d) 
(we use the same notations as Dubois and Lallemand [8]) 

mf ^^ = Gi {{Wj}o<j<d) ,d+\<i<N. (15) 



or 



mf ^^ = Gi {{Wj}^<j<d) = GtjWj = Gtjmj (16) 

In order to implement the linear stability analysis and recover the linearized macroscopic equations, we introduce 
the linearized form of Eq. (15) around reference states [3, 8, 6]. Using Eq. (16), the linearized description of Eqs. (9) 
and (10) can be written as follows 

mi(x, t + dt) = MiiMJ^"^ prMrix - ViSt, t\ (17) 

where the matrix ^ has the following form 



^ 






i,0<j<d 

{iij-Sij)^ 



(18) 



d+l<i<N,d+l<j<N 7 



where 0/^ = SfGij. A is matrix with the size (A/^ - J - 1) x (J + 1)), lis an identity matrix and the diagonal matrix S is 
defined by 

S = Diag(a^_0, Sd^u. ^ ..,SM^ . (19) 

d+l N-d-l 

In Eq. (17), the indices /, p and r indicate the summations from to N. 

In order to derive the linearized high-order equations, one assumes that the discrete single particle distribution fi 
belongs to C^iT x X) (a functional set, in which the element possesses a sufficiently smooth property with respect to 
the time domain T and spatial domain X ). This assumption is also used by Junk et. al for asymptotic analysis of the 
LBM [12]. This regularity hypothesis indicates that macroscopic quantities m/ are smooth ones and that the linearized 
system (17) is well defined. 

The next step consists of performing the Taylor series expansion of the right hand of Eq. (17), yielding 

mtix, t^St) = 2] —Mui-v^d^TMip'^^prmr, (20) 

where a indicates the summation from 1 to d. 

Now, we define the matrix A* = (A*. ) as follows 

A*r,n = ^Mui-v'^darMJ^'^,,. (21) 

When we need to derive equivalent equations or modified equations, it is difficult to use the matrix A*^ to carry out 
the calculations. In order to overcome this difficulty, we use the differential operators in spectral space. Let us note 
da = ika, with ka the wave-number in the a -direction and i^ = -1. Then, in spectral space, the matrix A*„ has the 
following form (A,„ = (A,7,„)o</<a^,o<;<a^) 

Air,n = ^Mui-iv^kaTMip'^pr- (22) 

Therefore, Eq. (20) can be rewritten as follows 

/-I 
m,(jc, t-\-St) = ^ dfAir^nrrir + Oidt-^). (23) 

«=0 



In order to derive the L-NSE corresponding to the L-MRT-LBM defined by Eq. (23), we introduce an original 
recursive algorithm. Given m/ = Wi(0 < i < d) (macroscopic conservative quantities), the algorithm is given as 
follows 

• Initial step. The initial O i and Bi are given as follows 

0,,;i = ^,,(0 < / < d\ 0,,;i = -^./J + 1 < / < TV), (24) 

Bij,l =Air,l^rj,l- (25) 

Let W = {Wi}o<i<d and m = {m/}o</<A^ denote the vector of the conservative quantities and the vector of all 
macroscopic quantities respectively. 

At the first order of St, for all macroscopic quantities, we have 

mi = O^i Wj + 0(St), 0<i<N. (26) 

By the matrix form, we have 

m = ^,iW -\- 0(St). (27) 

At the first-order of St, for conservative quantities, we have 

dtWi = Atr,l^rj,lWj + 0(St). (28) 

The matrix form is 

d,W = A iO,i Wy + 0(St). (29) 

• Recursive formula for all macroscopic quantities. 0^„ can be given as follows 

1 '^"^ St^ ^~^ 

^ij,n = -Ci'ij - Yj T^ir,n-lB[jn-l + ^ ^^^ ^ir,l^rj,n-l\ d+\<i<N (30) 

and 

^ij,n=Sij,(S)<i<d). (31) 

Eliminating the higher-order term of St^~^ , we have 

n-l 

^ij^n = Yj ^^^Coeff'(0,,;„ St, I), (32) 

where Coeff'(-, •, •) is a function which extracts the coefficients of the polynomials, for example, f(x) = 2/Lo ^i^^ 

CoQf^(f(x),x,i) = ai. (33) 

According to Eqs.(30) and (32), we have 

m = 0,^W-hO(^f). (34) 

• Recursive formula for conservative quantities. 5/^ „ is presented as follows 

^'> = - Z (7TI)! C-/ + Z St'-%„<P,j,„,,.,. (35) 

Eliminating the higher-order term of St'^'^, we have 

n-l 

Bij,n = Yj ^t^Coeff(Bij,n, St, /). (36) 

Now, for the conservative quantities, we have the following equation system 

dtW = B^n'W-^ OiSt""). (37) 

Using Eqs. (32), (34), (36) and (37), we can get the coefficient matrix of the conservative quantities at any order 
of St . Details are displayed in Appendix A. 



2.3. Illustrating example: application to a 2D MRT-LBM 

In this section, we illustrate the algorithm presented in 2.2 considering a 2D MRT-LBM. For the standard 2D 
MRT-LBM, the equilibrium distribution functions are described as follow [6, 13] 
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(38) 



where j;^ and j^y denote the x-momentum and y-momentum respectively, and p represents the density (Wo = mo = 
p,Wi = mi = jx,W2 = m2 = jy). The corresponding matrix ^ is given by 
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(39) 



The diagonal elements of the corresponding diagonal matrix S are set as follows 

'^O = '^l = ^^2 = 0, ^3 = Se, S4 = S^, 85 = 8^ = Sq, S^ = S^ = Sy. 



(40) 



For the original MRT-LBM [6], only Sy is a free parameter, and Se = 1.64, s^ = 1.54, Sq = 1.9. The analogous form of 
^ can be found in existing literature [6, 8, 13]. However, there still exist some slight differences. The derivation of ^ 
can be achieved by means of the first-order Taylor series expansion with respect to p, jx and jy at reference states. In 
the expression of ^, ([/, V) to denote the uniform flow velocity components. 
For the sake of convenience, we introduce the following relation 



1 



^^7= - 



(41) 



where rj stands for any notations in the set 



[e,€,q,y}. 



2.3.1. Considering the zero-mean flows (U, V) = (0, 0) 

Now, when the truncated error term is equal to 0(St^) , the coefficient matrix B 5 with the zero-mean flow can be 
described by the summation of five matrices. The first two matrices are given as follows, which describe the specific 
terms in the Navier-Stokes equations. 

The coefficient matrix associated with St^ is 






-k. 
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3 -^ 
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The coefficient matrix associated with St is given by 











U Q Kv Cr e Q '^x ^v Q '^v ^v 



(42) 







(43) 
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At the higher-order truncated errors of St, the coefficient matrices are given in Appendix B. Compared with 
the results given by Dubois and Lallemand [8], it is shown that the proposed algorithm yields the correct results. 
Meanwhile, the higher-order terms of St are off'ered by our algorithm explicitly. 



2.3.2. Considering the uniform flow (U,V) (U i^ or V i^ 0) 

When the truncated error term is equal to 0(6 f') , the coefficient matrix B2 with the uniform flow can be described 
by the summation of two matrices. The coefficient matrix associated with df is 



-\k^ + k^U^+kyUV -2k^U-kyV -kyU 

-\ky + kyV^ + k^UV -k^V -2kyV-k^U 



(44) 



The coefficients of St are given in Appendix C. 

From the coefficient matrix, it is clear that the correct convection terms of the L-NSE can be given by the L- 
MRT-LBM. However, the correct dissipation coefficients can not be obtained by the L-MRT-LBM with respect to 
the uniform flow, yielding the definition of a flow-dependent viscosity. Furthermore, this dependence also becomes a 
source of the non-Galilean invariance. 

From 2.3.1 and 2.3.2, it is known at the zeroth-order and the first-order of St, the relaxation parameters s^ and Sq 
have no influence on the recovered L-NSE. Here, the given form of the L-NSE is generated with respect to the small 
perturbations of the density p and the momentum quantities (jx, jy). 

3. Optimization strategies of free parameters in the MRT-LBM 

In this section, the original optimization strategies of free parameters in the MRT-LBM are proposed based on 
the matrix perturbation theory and the modified equations. The optimized parameters will be determined in order to 
obtain the optimal dispersion/dissipation relations. 

3.1. The matrix perturbation theory for the L-NSE corresponding to the L-MRT-LBM 

From Eq. (37), it is known that the dispersion and dissipation relations of the L-MRT-LBM are determined by the 
matrix En e c^'^+^M^+d)^ where c'^^+dM^+d) denotes the (1 -h J) x (1 -h d) complex matrix set, if the truncation error is 
developed up to the ^th-order of St. Here, B refers to the coefficient matrix of the exact L-NSE in wave-number space 
similar to the matrix En in Eq. (37). This means that for the exact L-NSE, one has the following expression 

dtW = B' W. (45) 

For a given n, if the errors terms with order higher that St^ are neglected, the main deviation of dispersion and 
dissipation relations between the L-MRT-LBM and the L-NSE originates in the diff'erences between eigenvalues of 
Bn and those of B. 

Now, we introduce the perturbation matrix Mg e c^^+dM^+d) (defined as 

B^n = B-^ Ms. (46) 

Let B have eigenvalues Ai,. . .,An and Bn have eigenvalues li, . . . , 1„. The spectral variation of Bn with respect to B 
is [14] 

SYB(B,n) = max min \Ai - Aj\. (47) 

i J 

Then [14], 

SYB(B,n) < (\\B\\2 + \\Bj2)'-'^^'^'^\\Ms\\l^^'^'\ (48) 

where || • II2 is the spectral norm of matrices. For all A e c^d+i)x(d+i) ^ y . y^ ^^ defined by 



IIAII2 = V^max(A^A) = cr^ax(A), (49) 

where A^ denotes the conjugate transpose of A. Aj^axiA^A) is the largest eigenvalue of A^A and (Tniax(A) is the largest 
singular value of A. In order to establish the direct relation between elements and eigenvalues of matrices explicitly, 
the Frobenius norm is given as follows for any A = (<3/y)o</<j,o<j<7 ^ c^d+i)x(d+i) ^^j^ 



\\A\\f = 



d 
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A Z Z l^^-^l' = Vtrace(A^A) = a 2, ^?' (^0) 



where ct/ denotes the singular values of A. 

Furthermore, let cr(B) = {Aq, . . . , Ad+i}, the multiset of B's eigenvalues, and set 

A = diag(/lo, . . . , Ad+i), Ar = dmg(Ar(o), . . . , Ar(d+i)), (51) 

where r is a permutation of {1, . . . , J + 1}. Let cr(Bn) = {Aq, . . . , A^+i}, the multiset of En's eigenvalues, and set 

A = diag(lo, . . . , ^d+i), K = dmg(Ar(0), . . . , Ar(d+i)), (52) 

Then, there exists a permutation r such that the following inequality is satisfied [14, 16] 



I|A-Arll2<2 



J+1 



(ll^lb + \\BJ2)'-'^^'^'^\M,\\\^^'^'\ (53) 



Since ||Me||2 < ||Me||/7, the minimization of SVfi(5,„) or ||A - A-^lb means ||Me||/7 should be minimized to reduce both 
dispersion and dissipation errors associated with the MRT-LBM schemes. 

3.2. Optimization methodology 

The following wave-number definition in Eq.(37) is considered 

k^ = k- cos(^), ky = k- sm(6). (54) 

Substituting Eq. (54) into Eq. (37) and considering the uniform flow (U, V), we can get the following formal expres- 
sion for Bn 

n-\ 
B,n = Yj ^t'k'^'Hj,n(^^ S, U, y), (55) 

standard MRT-LBM [13]) 



where b^.. is a function of 6, S and (U,V). S denotes the following set (about extra parameters, refer to the non- 



S = {a,l3. A, (Te, (Te, (T^, CFy}. (56) 

For the standard MRT-LBM, the parameters a, ft and A are equal to 1, -3 and -2, respectively. In order to handle 
the influences of the uniform flows, we consider the following uniform flows 

U = Urn' cos(?^), V = Urn' sin(?^), (57) 

where Um denotes the magnitude of the uniform velocity. Now, Eq. (55) has the following form 

n-l 

B,, = Y,^t^k^^'b\jn(e,i^,E,u^). (58) 

Furthermore, Eq. (58) is rewritten as follows 

1 1 '^ 

Bij,n = -^ijA^tK 0, ^, S, u^) = - Yj^6tkib\r:l^{e, ^, S, uj. (59) 

For the L-NSE, there exists the similar matrix ^ defined by ^ = 6tB. It is known that for the MRT-LBM, the 
dispersion error corresponding to the L-MRT-LBM stems from the odd-order coefficient matrix Mf of dtk in Me = 
StMg = B^n - ^, and the dissipation error comes from the even-order coefficient matrix Mf of dtk in M^. So, the 
perturbation matrix M^ can be expressed as follows 

Ate = M? -h Mf . (60) 



3.2.1. The zero -mean flow case 

Considering Um = 0, then, Me = ^,„ - ^ is a function of Stk, 6 and H. According to Sec. 2.3.1, 

n n 

Me = Yji^tk)'b'i-]l{0, &, E, u„) = Y,m)'b\-l{0, H), (61) 

1=3 1=3 

Mf= Yj (.Stkf^'bl^{e,E), (62) 

\<l<n,2l+\<n 

Mf= Yj (.6tkfbfr^\0,E). (63) 

\<l<n,2l<n 

According to the theory of the finite diff'erence method (FDM) [17], Eq. (37) can be regarded as the modified 
equation of Eq. (45). In modified equations, the higher even-order derivatives beyond Eq. (45) cause numerical 
dissipation and the higher odd-order derivatives cause numerical dispersion [17]. 

In order to reduce the dispersion error, when cr^ and CFy are specified, it is proposed here to minimize the following 
cost function: 

F%K)= r { \\M%l&e&{dtk). (64) 

Jo Jo 

For the standard MRT-LBM, the parameters cr^ and CTq need to be determined. The corresponding conditions are 

cr, >0,cr, >0, (65) 

that is to say, 

^e,^^e(0,2]. (66) 

In order to separate the kinetic modes form the modes directly aff'ecting hydrodynamic transport, Lallemand and Luo 
[6] suggested that s^ and Sq should be kept slightly larger than 1. Accordingly, it was implied that Se,Sq e (1,2). In 
this paper, s^, Sq are taken in the range (0, 2]. 

Furthermore, the same method can be used to reduce dissipation error. The corresponding cost function is 

F\E) = r ['' IIMf ll^d^d(M). (67) 

Jo Jo 

If both of dispersion and dissipation errors need to be reduced, according to Eqs. (48) and (53), minimizing error 
between hydrodynamic modes of the L-NSE and the L-MRT-LBM is achieved by minimizing 

r(H) = r ['' \\Me\\ldOd(Stk) = r f \\\mX + \\M^\\l)ded(Stk). (68) 

Jo Jo Jo Jo 

3.2.2. The non-zero mean flow case 

Considering Um i^ 0, then. Me = ^,„ - ^ is a function of dtk, 6, ?^, S and Um- It is known that it is difficult for the 
non-zero mean flows to determine the values of free parameters locally, because the optimization problems must be 
solved at each lattice node. In order to avoid solving optimization problems locally, the minimization problems will 
be integrated with respect to ?^ and Um- Here, it is necessary to mention that the similar relations of Eqs. (62) and (63) 
about Mf and Mf are not satisfied for the non-zero mean flows. 

In order to minimize the dispersion error, the following cost function is introduced 

G^(S)= ||Mf|||d^d(M)dMw^, (69) 

Jo Jo Jo Jo 

where wq is the upper bound of lattice velocity magnitude. Generally, u^ is taken equal to 0.2. Similarly, in order to 
minimize the dissipation error, we have the cost function 

G\E)= ||Mf||2dM(5ffc)dMM„. (70) 

Jo Jo Jo Jo 



For the non-zero mean flows, the optimization problem of minimizing error between hydrodynamic modes of the 
L-NSE and the L-MRT-LBM is associated with the foUowing cost function 



r*UQ rO.n r*7: pin 

^(S)= \\M£pdeA{6tk)AMu. 

Jo Jo Jo Jo 



(71) 



3.2.3. A non-zero mean flow case and separating the bulk-viscosity terms from the dissipation coefficient matrix 

It is observed that when the shear viscosity is very smaU, the magnitude of the bulk viscosity is very sensitive 
to the stability numerically and theoretically. If the bulk viscosity is also too small, the MRT-LBM schemes will be 
unstable. Meanwhile, although we adopt the optimization strategy detailed in Sec. 3.2.2 and the optimized MRT-LBM 
appears to be more stable than the original MRT-LBM, the stability of the obtained MRT-LBM is still very sensitive. 
Furthermore, it is known that for linear acoustic problems, the values of shear and bulk viscosity are often very 
small and the dissipation efl'ects from the shear and bulk viscosity can nearly be neglected. In the simulations, if the 
bulk viscosity is too large, the pressure fluctuations will be damped very significantly. This over-dissipative behavior 
should be avoided for aeroacoustic problems. In order to handle the low bulk viscosity problems, we propose a new 
optimization strategy. For the uniform flows, the linearized convection terms are given by the matrix (44) in spectral 
space, and the recovered linearized dissipation matrix is also given by a matrix in Appendix C. The exact linearized 
dissipation coefficient matrix with the uniform flows is given by (the simple derivations with respect to the perturbation 
of p, jx and jy are neglected) 











\J o /Cj^ Cr £ o K/X CTy o Jv-y CTy 







O l\'X '*'V Cf £ 





^ K'X i^y ^ e 



(72) 















-\k^cre 


— l^Kxtiy Cq 





— l^KxI^y C £ 


-\ky^(Te _ 



It is clear that the coefficient matrix (72) has the same expression as the coefficient matrix (43). The coefficient matrix 
corresponding to bulk viscosity is given by 



(73) 



The new perturbation matrix Me is defined by Me - ^,n - ^ + ^. The matrix Me possesses the information of 
the bulk viscosity at the first-order of dt. The optimization strategies are kept the same as those in Sec. 3.2.2. The 
parameter (T^, which is related to the bulk viscosity, can be taken as an user- specified or free parameter for acoustic 
problems. If the MRT-LBM is considered as a high-precision solver for the nearly incompressible flows, (Je can be 
set as a free parameter. 

4. von Neumann analysis of the dispersion and dissipation relations 

In this part, the optimization strategies will be investigated by means of von Neumann analysis. At the same time, 
we will give some numerical results of the dispersion and dissipation relations. 



4.1. Theoretical dispersion and dissipation relations for the L-MRT-LBM and the L-NSE 

In order to apply von Neumann analysis to validate the optimized parameters, it is necessary to give the expressions 
of the L-MRT-LBM in frequency- wave number space. Considering a uniform mean part f^ and a fluctuating part f/, 
the equilibrium distribution function can be linearized as [3, 4] 



f:'^\{ff+f/h<_j,N)=fr'"+ 



df 



(eq) 



dfj 



fj=ff 



■fj' + 0{f'^). 



(74) 
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Table 1 : Optimized free parameters for zero-mean flow case. 



Groups Methods Un 



F%E) 



r(E) 



A Min. (64) 0.0025 

B Min. (68) 0.0025 

C Min. (64) 0.1 

D Min. (68) 0.1 



105.468091254867 

105.465307838135 

26.3631592758091 

26.3520430827600 



17.9024342612509066 

\ 

17.9107477965877778 

\ 



\ 
17.9030645832220686 

\ 
17.9208148202264042 



Then, considering a plane wave solution of the linearized equation 

fj = AjQxp[i(k ' X - ojt)], 



(75) 



according to Eqs. (74), (75) and (14), we get the following eigenvalue problem for the L-MRT-LBM in frequency 
wave number space [3, 4, 6] 



where the matrix Af^ = A'^ [/ - M-^SMN^^^], and A^^^^ is defined by 



A^'.^^ 



fj 



(76) 



(77) 



fj-^ 



(78) 



For the L-NSE, the analytical acoustic modes oj- and shear modes oj^k) are given by [18] 

^ Re[oj^(k)] = |k|(±c,(k) + |u|cos(k^)), 
Im[oj^(k)] = -|k|2i (^y(k) + rj(k))) , 

Re[oj'(k)] = |k||u|cos(k • u), 
Im[oj'(k)] = -|kpy(k), 

where v is the shear viscosity and t] is the bulk viscosity. 

4.2. Optimized free parameters with the zero-mean flows or with the non-zero mean flows 

First, we set ^ = 5 and t/^^ = in Eqs. (37) and (59). When cTe = 0.0025 and o-y = 0.0025, the analytic expressions 
of the problems (64) and (67) are given in Appendix D. The optimized results are given in Table 1. 

It is observed that f or ^ = 5 and Um = 0, the optimized value of cr^ is close to zero or equal to zero, this is to say, 
the corresponding s^ is close or equal to 2. When Ce and o-y are attenuated, numerically, cr^ is equal to 0. However, 
the magnitude of cr^ increases. Furthermore, the values of F^(E) and 'F(S) are around 17.9. In Fig. 1, the 3D surfaces 
of F^(E) are displayed. From the figure, it is known that the function F^(E) is convex and the extreme value can be 
found along a^ = 0. According to present numerical investigations, the value of a^ is always close or equal to 0. For 
the zero-mean flows, when o-g and o-q are specified, we can specify o-^ equal to or slightly larger than in order to 
simplify the minimization problems for practical applications. 

The Fig. 2 displays the dispersion and dissipation profiles of the L-MRT-LBM. In the original MRT-LBM, the 
free parameters cr^ and o-q are taken equal to the recommended values [6] and the optimized free parameters are given 
by Group A in Table 1. From numerical results in Fig. 2, it is observed that when the wavenumber k is perpendicular 
or parallel to x-axis, the optimized MRT-LBM has the same dispersion relations as the original MRT-LBM. When 
6 is equal to other values, the optimized MRT-LBM performs better than the original MRT-LBM. These results also 
indicate that the cross derivatives in Eq. (37) are the main source of dispersion error. For the dissipation relations, 
it is observed that by the minimization problem, we can enhance the stability of the MRT-LBM. From the profiles 
of the dissipation relations, there exists one unstable mode the imaginary part Im(^) of which is larger than 0, when 
0.75 < kAx < 1.75, e^OsindO^ 7t/2. 

According to the definition of Me in Sec. 3.2.2, we now consider n = 4 and Um = 0.1 in Eqs. (37) and (59). Under 
these conditions, the dispersion and dissipation relations are investigated firstly. In Fig. (3), the 3D surfaces of ^(H) 
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Figure 1: The 3D surfaces of the minimization function F^(E): (left) o-g = 0.0025, cry = 0.0025 (right) (Te = 0.1, cry = 0.1 
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are shown. It is discovered that the function Q(E) is convex. Some optimized results for specified o-g and CTy are given 
in Table 2. 

In Fig.4, we show the dispersion and dissipation profiles. It is seen that the best shear mode description is given 
by the optimization problem ^(S). At the same time, the optimized MRT-LBM is more stable than the original MRT- 
LBM. However, it is discovered that the dispersion error is not improved for ^ = 4 by the minimization problem 
(71). If we want to reduce the influence of the non-zero mean flow on the dissipation relation and avoid handling 
lengthy mathematical expressions, it is suitable to choose ^ = 4 in Eq. (37). Furthermore, when = n, the optimized 
MRT-LBM has the similar dispersion and dissipation profiles with the original MRT-LBM in Fig. (4)-b. According to 
authors' numerical investigations, when the angle 6 is equal to 0, 7r/2 and n, there always exist the similar dispersion 
and dissipation profiles between the optimized MRT-LBM and the original MRT-LBM. Based on the definition of Ale 
in Sec. 3.2.2, the results of numerical studies have shown that by the optimization problems (69), (70) and (71), the 
dispersion error was not reduced when ^ < 5 in Eq. (37) for the uniform flows. 

According to the definition of Al^ in Sec. 3.2.3, we consider n = 4 and Um = OA in Eqs. (37) and (59). For very 
small shear and bulk viscosity parameters, we show some results in Table 3. 

In Figs. 5 and 6, the dispersion and dissipation relations are given for both the optimized MRT-LBM and the 
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Figure 2: Dispersion and dissipation profiles of the L-MRT-LBM based on recommended free parameters (except Se and Sy) [6] and optimized 
free parameters (group- A in Table 1). The relaxation parameters Se and Sy are kept equal to those of the original MRT-LBM and the optimized 
MRT-LBM. The (#-1) figures display the dispersion profiles (magnified locally) and the straight lines represent the exact dispersion solutions. 
The (#-2) figures display the dissipation curves and the line is the expected shear mode and acoustic mode dissipation. The angle 6 between the 
wavenumber k and the x-axis: (a) 6 = 0;(b) 6 = n/6; (c) 6 = n/4; (d) 6 = n/2; (e) 6 = 2n/3; (f) 6 = 3n/4. (The symbol "#" stands for the characters 
a, b, c, d, e and f.) 



Table 2: Optimized free parameters for uniform flows (u,n = 0.1) based on the perturbation matrix Me in Sec. 3.2.2. The parameters are obtained 
by the minimization (71) 



Groups 



(Te = (Ty 



^(S) 



A 0.001 
B 0.0025 



0.00751873323089156 
0.0187982349323006 



0.00171909400064198 
0.00429903714246050 



14.4537555316616474 
14.4550408968152340 
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Figure 3: The 3D surfaces of the minimization function T^(S): (left) o^e = 0.001, o-y = 0.001; (right) cTe = 0.0025, o-y = 0.0025. 
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Figure 4: Dispersion and dissipation profiles of the L-MRT-LBM based on recommended free parameters (except Se and Sy) [6] and optimized 
free parameters (group- A in Table 2). The relaxation parameters Se and Sy are kept the same values for the original MRT-LBM and the optimized 
MRT-LBM. The (#-1) figures indicate the dispersion profiles (magnified locally) and the straight lines represent the exact dispersion solutions. The 
(#-2) figures indicate the dissipation curves and the line is the expected shear mode and acoustic mode dissipation. The angle ?^ = n/4, U = O.l and 
V = 0.1. The angle k • u between the wavenumber k and u: (a) k • u = 0, ^ = 7r/4;(b) k • u = n/4, 6 = n; (c) k • u = n/2, 6 = 37r/4. (The symbol 
"#" stands for the characters a, b, and c.) 



Table 3: Optimized free parameters for uniform flows (um 
by the minimization (71). 



: 0.1) based on the perturbation matrix Me in Sec. 3.2.3. The parameters are obtained 



Group 


O-e 


(Ty 


o-e 


O-q 


^(S) 


A 
B 


0.0025125628 
0.0000025 


0.00001 
0.00001 


0.00947095595580122 
0.0000311837523990053 


0.00181622011980382 
0.00000760552251906507 


14.453655614637361 
14.45351071030945 



original MRT-LBM. It is discovered that, by the new definition of the perturbation matrix in Sec. 3.2.3, the optimized 
MRT-LBM is always stable for the very small viscosity and bulk viscosity. Furthermore, when the bulk viscosity 
(corresponding to Fig. 6) is very small, the numerical shear modes given by the optimized MRT-LBM agree with the 
exact shear modes very well. The obtained shear modes are nearly exact compared with the shear modes given by the 
original MRT-LBM. In order to observe the details of the optimized dissipation relations and the exact relations, in Fig. 
7, the locally-magnified dissipation relations are given. It is clear that we observe a lower dissipation of the acoustic 
modes for the optimized MRT-LBM. From Figs. 5 and 6, it is also clear that when the bulk viscosity become smaller, 
the modes of the original MRT-LBM become more unstable for all values of the angle 0. These results indicate that 
the original MRT-LBM is not suitable for aeroacoustic problems, because the bulk viscosity in the original MRT-LBM 
can not be chosen to be arbitrarily small. This limitation in the original MRT-LBM means that there exists a strong 
dissipation of the acoustic waves in the numerical simulations. The new definition of the perturbation matrix Me in 
Sec. 3.2.3 coupled with the optimization strategy (71) overcomes this drawback of the MRT-LBM under the premise 
of guaranteeing the stability. 

5. Numerical simulations of acoustic problems 

In this section, the classical acoustic problems will be simulated by optimized MRT-LBM. At the same time, some 
comparisons between the D2RP MRT-LBM and the original MRT-LBM are given. 
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Figure 5: Dispersion and dissipation profiles of the L-MRT-LBM based on recommended free parameters (except Se and Sy) [6] and optimized 
free parameters (group- A in Table 3). The relaxation parameters Se and Sy are kept the same values for the original MRT-LBM and the optimized 
MRT-LBM. The (#-1) figures indicate the dispersion profiles (magnified locally) and the straight lines represent the exact dispersion solutions. The 
(#-2) figures indicate the dissipation curves and the lines are the expected shear mode and acoustic mode dissipation. The angle d^ = n/4, U = 0.1 
and V = 0.0. The angle k • u between the wavenumber k and u: (a) k • u = 0, ^ = 0;(b) k • u = n/3, 6 = n/3; (c) k • u = n/4, 6 = n/4; (d) 
k • u = n/2, 6 = n/2. (The symbol "#" stands for the characters a, b, c and d.) 





Table 4: Specified relaxation parameters and optimized free relaxation parameters for zero-mean flows. 




Groups 


Se = Sy S^ Sq V p 


Ps 


A 
B 


1.99044751 2 0.00875438872 0.000799861 1 
1.95321 2 0.04126919093 0.00399257 1 


0.01 
0.01 



5.1. Acoustic point source 

In this part, we validate the optimized MRT-LBM by an acoustic point which sends out a sinusoidal signal [19, 20]. 
The point source is set by the following density configuration [19] 



p(x,0 =po+P.sin(y^|, 



(79) 



where ps is the point source amplitude, and T the period of the oscillation with respect to lattice units. In order to 
avoid nonlinear wave eff'ects, it is necessary thatp^ «: po. The macroscopic velocity (w, v) at the point source is equal 
toO. 

It is known that the sound speed of the D2Q9 MRT-LBM is equal to c^ = 1/ V3 = 0.57735. In order to avoid 
the eff'ects of boundaries, the wave propagation is limited in the computational domain D = 100 x 100. The lattice 
nodes are 101 x 101 and the acoustic point source is set in the center of computational domain. So, if the lattice 
computational time t is in the range (0, L50/c^J] = (0, 86], the wave will be limited in the domain Q. Now, we will use 
parameters given in Table 4 to validate the optimized MRT-LBM with the periodic boundary conditions. 



Table 5: Specified parameters and optimized free 
parameters are obtained by the minimization (71) 


parameters for uniform flows based on tlie perturbation matrix Me 


in Sec. 


3.2.3. Tlie relaxation 


Groups 


Se 


Sy 


Se 




^« 


A 
B 


1.99 
1.99999 


1.999960001 
1.999960001 


1.962820428 
1.999875273 




1.992761413 
1.999969578 
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Figure 6: Dispersion and dissipation profiles of the L-MRT-LBM based on recommended free parameters (except Se and Sy) [6] and optimized 
free parameters (group-B in Table 3). The relaxation parameters Se and Sy are kept the same values for the original MRT-LBM and the optimized 
MRT-LBM. The (#-1) figures indicate the dispersion profiles (magnified locally) and the straight lines represent the exact dispersion solutions. The 
(#-2) figures indicate the dissipation curves and the lines are the expected shear mode and acoustic mode dissipation. The angle ?^ = n/4, U = 0.1 
and V = 0.0. The angle k • u between the wavenumber k and u: (a) k • u = 0, ^ = 0;(b) k • u = n/3, 6 = n/3; (c) k • u = n/4, 6 = n/4; (d) 
k • u = n/2, 6 = n/2. (The symbol "#" stands for the characters a, b, c and d.) 



In Figs. 8 and 9, the density contours and 3D surfaces are shown. Figs. 8(a) and 9(a) show that the optimized 
MRT-LBM corrects the anisotropy and annihilate the spurious fluctuations of waves. These results are better than 
the BGK-LBM and the original MRT-LBM. In Figs. 10 and 11, the waves along the line y = 51 are given for three 
diflferent methods at ^ = 80 and t = 100. It is clear that the optimized MRT-LBM is more eflfective than the BGK-LBM 
and the original MRT-LBM. 

The numerical results demonstrate that by determining the free parameters, we can reduce the dispersion error and 
the isotropy error of the MRT-LBM. 

5.2. Acoustic pressure pulse 

In this part, the quality of the optimized MRT-LBM is assessed considering the acoustic pressure pulse problem. 
These comparisons provide an evidence for the accuracy of computed solutions. In order to compare the results with 
the exact solution and neglect the efl'ects of the dissipation, the very small values of shear and bulk viscosity are chosen 
for a MRT-LBM. The initial perturbation is given by a Gaussian density distribution at the center of the domain at 
t = 0[n] 

( p(x, y, 0) = po + sexp(-ar]^), 

I u(x,y,0) = Uo, (80) 

[ v(x,j,0) =0, 

where a is related to the half- width Gaussian , b, by a = ln2/b^. t] is defined by 7/ = \(x- Uotf- + j^ , which is 
equal to the radial coordinate at ^ = 0. e is the density pulse amplitude. The analytical solution of the problem for the 
pressure and density can be given by a zero-order Bessel function /q [1 1] 



p(x,yj) =po + 



2a 



f 

Jo 



cxp(-e/(4a))cos(Cs^t)Jo(^ri)^d^. 



(81) 



It is noted that Eq. (81) does not include the dissipation introduced by the viscosity of the fluid. So, in order to 
implement the computation, the influence of viscosity on pressure wave must be minimized by the small magnitude of 
viscosity. The computational domain is [0, 1] x [0, 1]. The parameters s and b are equal to 0.01 and 0.04, respectively. 
In Figs. 12 and 13, we show the horizontal density profiles at y=0.5. The simulation physical time t = 0.4. From these 
figures, it is clear that the results by the optimized MRT-LBM are superior to the results by the original MRT-LBM. 
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Figure 7: Locally-magnified dissipation profiles: Dissipation profiles of exact solutions (acoustic mode: w-'^; shear mode: w^'^) and optimized 
free parameters (group-B in Table 3; acoustic mode: oj-; shear mode: oj^). The angle ?^ = n/4, U = O.l and V = 0.0. The angle k • u between the 
wavenumber k and u: (a) iTu = 0,6 = 0;(b) iTu = n/3, e = n/3; (c) iTu = n/4, 6 = n/4; (d) iTu = n/2, 6 = nil. 
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(a-1) t = 80 








(a-2) r = 80 




(b-2) t = 80 




(c-1) t = 80 



(c-2) r = 80 



Figure 8: The density contours (left figures) and the 3D surfaces (right figures) by the parameters of group-A in Table 4: (a) the BGK-LBM; (B) 
the original MRT-LBM; (C) the optimized MRT-LBM. 
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(a-2) t = 80 




(b-1) t = 80 



(b-2) t = 80 





(b-1) t = 80 



(b-2) t = 80 



Figure 9: The density contours (left figures) and the 3D surfaces (right figures) by the parameters of group-B in Table 4: (a) the BGK-LBM; (B) 
the original MRT-LBM; (C) the optimized MRT-LBM. 
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Figure 10: The comparisons of cross-profiles of the density p at y=51 (the lattice time t = 80): (a) Results obtained by the parameters of group-A 
in Table 4;(b) Results obtained by the parameters of group-B in Table 4. 
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Figure 11: The comparisons of cross-profiles of the density p at y=51 (the lattice time t = 100):(a) Results obtained by the parameters of group-A 
in Table 4;(b) Results obtained by the parameters of group-B in Table 4. 
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(a) Horizontal density profiles at y=0.5 



(b) Horizontal density profiles at y=0.5 



Figure 12: The comparisons of the cross profiles of the density p at y=0.5: The horizontal velocity Uq = 0. (a) The results of the optimized 
MRT-LBM are obtained by the parameters of group -A in Table 5 (b) The results of the optimized MRT-LBM are obtained by the parameters of 
group-B in Table 5. (The computational physical time t = 0.4) 



Obviously, the amplitudes of crests and troughs are damped by the bulk viscosity in the original MRT-LBM. It is 
demonstrated that by the proposed strategy (71) and the new definition of Me in Sec. 3.2.3, the dissipation influence 
from the bulk viscosity on pressure waves can be reduced to a negligible level. Meanwhile, the bulk viscosity can 
be attenuated. In Fig. 14, the figures of density distribution are given. In order to test the accuracy, the following 
L^-norm relative error is defined 



^L2 



A 



■^A^nodes /^th 



num\2 



i:;iV"(pf-pr"'> 

Y Anodes /^th\ 2 



■^/=1 



(Pfy 



(82) 



where pf and p^^^ denote the theoretical and numerical solutions, respectively. In Fig. 15, the L^-norm relative 
errors of density are given based on lattice nodes in a log-log coordinate. From this figure, it is discovered that the 
convergence orders of pressure pulse for the optimized MRT-LBM and the original MRT-LBM are very close to 1. 
However, the result given by the optimized MRT-LBM are better than that given by the original MRT-LBM. 

6. Conclusion 

In this paper, we have proposed several numerical strategies to reduce the dispersion/dissipation errors (regarded 
the optimized MRT-LBM as the D2PR-LBM). We also gave an easy-to-use algorithm to derive linearized Navier- 
Stokes with the high-order truncation errors starting from the linearized MRT-LBM. Von Neumann analysis of the 
isothermal linearized MRT-LBM and the linearized BGK-LBM has been made to investigate the dispersion and dis- 
sipation relations. The von Neumann stability analysis also shows that by optimizing free parameters, the acoustic 
modes are decoupled from the shear modes and other modes with respect to zero-mean flows when the truncation 
error is up to 0(St^). For uniform flows, it is discovered that when the truncation error is up to 0(St^), by optimization 
strategies, the dissipation error can only be reduced and the influences from mean flows on dissipation relations are 
also reduced. The stability of the MRT-LBM is enhanced. Especially, when the shear viscosity and bulk viscosity are 
very small, the optimized dissipation relation is nearly exact. The optimized MRT-LBM can annihilate the spurious 
waves and the isotropic error sufl'ered by the original MRT-LBM and reduce the over-damping influence of the bulk 
viscosity on pressure waves. Numerical simulations of acoustic problems demonstrated that for acoustic problems, 
the optimized MRT-LBM is more efl'ective than the original MRT-LBM. 



25 



1.0020 



1.0015 



1.0005 



1.0000 



0.9995 



0.9990 



- Exact 
Optimized MRT-LBM 
Original MRT-LBIVI 




1.0020 



1.0015 



1.0010 



1.0005 



1.0000 



0.9995 



0.9990 



- Exact 
Optimized MRT-LBM 
Original MRT-LBM 




(a) 



Horizontal density profiles at y=0.5 



(b) Horizontal density profiles at y=0.5 



Figure 13: The comparisons of the cross profiles of the density p at y=0.5: The horizontal velocity Uq = 0.1. (a) The results of the optimized 
MRT-LBM are obtained by the parameters of group -A in Table 5 (b) The results of the optimized MRT-LBM are obtained by the parameters of 
group-B in Table 5. (The computational physical time t = 0.4) 
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(a) Density distribution by the optimized MRT-LBM 



(b) Density profiles by the original MRT-LBM. 



Figure 14: The comparisons of density distribution with the horizontal velocity Uq = OA. (a) The results of the optimized MRT-LBM are obtained 
by the parameters of group-A in Table 5 (b) The results of the original MRT-LBM. (The computational physical time t = 0.4) 
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Figure 15: The comparisons of relative errors of the density p in the log-log coordinate. The results of the optimized MRT-LBM are obtained by 
the parameters of group-A in Table 5 (The computational physical time t = 0.4) 
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Appendix A. The Taylor expansion strategy of the L-MRT-LBM in wave-number spaces 

In this part, we show the derivation details from the L-MRT-LBM to the L-NSE and the truncation error is up to 

(1) When / = 1 in Eq. (23), we have 

Wi = m„ < / < J, 



mi = -^ijWj -h 0(St), d<i<N. 

Si 

We rewrite Eqs. (A.l) and (A.2) as a uniform expression 

m = ^ij,iWj -h 0(St), 0<i<N, 

where 

Om = Sij(0 < i < d\ O^i = -^ij{d +l<i<N). 



(2) When 7 = 2 in Eq. (23), we have 

When < / < fi?, we have 
Because A,>,o = <5,>, we have 
That is. 



nii + dtdtifij = Airfiifir + dtAjr^inir + 0(6r). 
Wi + 6td,Wi = Airflnir + StAir,imr + 0{6f-). 

Wi + 6td,Wi = SirMr + StAir,imr + 0{6t^). 

d,Wi = Air,imr + 0(6t). 
27 



(A.1) 

(A.2) 

(A.3) 
(A.4) 

(A.5) 

(A.6) 

(A.7) 
(A.8) 



According to Eq. (A. 3), we have 

dtWi = Atr,l^rj,lWj + 0(St). (A.9) 

Let 

Bij,l = Air,l^rj,l- (A.IO) 

Then, we obtain 

dtWi = Bij^iWj ^ 0(St). (A.ll) 

When d -\- 1 < i < N,WQ have 

rui + Stdtrui = A.y^o^r + StAir^irUr + 0(St^). (A. 12) 

It is known that when d -\- I < i < N, A^y^o is defined by [8] 

Air,omr = (Str " ^.v)^, - "i^ijWj. (A.13) 

So, we have (the combination of Eqs. (A.13) and (A.2)) 

nii + St^ij^idtWj = (Sir - Sir)mr - ^ijWj + dtAir,l^rj,lWj + 0(^^^), (A.14) 

nii = -(-^.7 - ^^0,,,i%i + StAtr,i^rj,i)Wj + 0(^r2). (A.15) 

Now, introducing 0/y;2 as follows 

^ij,2 = ^ij(0 < i < d\ (A. 16) 

we have 

0.7,2 = -(-Smk,iBkj,i + ^.7 + 6tAir,i^rj,i), (^ + 1 < / < A^). (A.17) 



mt = ^.2Wj + 0(St^). (A.18) 



So, we obtain 

(3) When / = 3 in Eq. (23), we have 

Sf o n 

Mi + ^^^^m/ + -—dtnii = Airfiifir + ^M/^jm^ + ^rA/^^2^r + 0(^r). (A. 19) 

When < / < J, we have 

dtWi + ^d,Wt = Air^iMr + dtAir,2mr + 0(^^2). (A.20) 

By Eqs. (A.ll), (A.3) and (A.18), we get 

dtWi = -f^Btj,lWj + Atr,l^rj,2Wj + StAtr,2^rj,lWj + OiSt^ (A.21) 



Let 

^.7,2 = -^Bij^iWj + A,v,lO,y,2W; + StAir,2^rj,U (A.22) 

we get 

^.^/ = Btj^2Wj + O(^r^) (A.23) 

When J + 1 < / < A^, by Eqs.(A.18), (A.3) and (A.23), we have 

1 Sf o Sf . 

Mi = -{-dt^ik,2Bkj,2 - —^ik,lBlji + ^,7 + StAir,l^rj,2 + -^Atr,2^rj,l)Wj + 0(St^). (A.24) 

Now, introducing 0.73, we get 

0,7,3 = ^.7(0 < / < d), 

28 



and for ( J + 1 < / < A^) 

1 Sfi . o 

^ij,3 = -{-^t^iKlBkj,! - -^^ik,lBlj^ + ^,7 + StAir,l^rj,2 + StAir,2^rj,l)- (A.25) 

Si Z ! •^' 

In order to restrict the truncated error of Eq. (A.25) equal to 0(Sp), we rewrite Eq. (A.25) as follows 

2 



/=0 

So, we have 

(4) When / = 4 in Eq. (23), we have 



07,3 = ^^^^Coeff(0^3,^^/). 

mi = ^ij,3Wj + 0(St^). (A.26) 



Mi + Stdttfii + -—d^nii + -—d^nii = Air^nir + ^M/^jm^ + dt Air^2mr + ^rA^y^m^ + 0{dt ). {A21) 

When < / < J, we have 

dt^i + ^^?W/ + -^d]Wi = Air^iMr + StAir^Mr + ^^^A,r,3m, + O(^r^). (A.28) 

By Eqs. (A.26), (A.23) and (A. 11), we get the R.H.S of Eq. (A.28), 

R.H.S = Air,i^rj,3Wj + dtAir,2^rj,2'^j + ^^^A.^B^r;,! W,- + Oidt"). (A.29) 

By Eqs. (A.23) and (A. 11), we get the L.H.S of Eq. (A.28), 

L.H.S = d,Wi + l^^^,- + f^4i W,-. (A.30) 



So, we have 

dtWi = (-^Bf.^ - ^—B',., + Air,l^rj,3 + ^M,v,20,,;2 + St'Air,3^rj,l)Wj + 0(^^^). (A.31) 

Let 

St . Sf . o 

^aS = -^^^2 - ^^ai + ^^>,l^rj;3 + StAir,2^rj,2 + ^rA,v,30,y,i, (A.32) 

we restrict the truncated error of Eq. (A.32) equal to 0(Sp) and get 

2 



St o St^ . 



/=0 

Then, we have 



When J + 1 < / < A^, we have 

St^ ^o St^ 



Bij3 = ^ ^^^Coeff (5,7,3, St, /). (A.33) 

dtWi = Bij^sWj ^ 0(St^). (A.34) 



rrii + StdtMi + ^dfrrii + ^dftfii = Air^mr + StAir^irUr + St^Air^2mr + St^Air^^^r + 0(^^^). (A.35) 

By Eqs. (A.3), (A.ll), (A.18), (A.23), (A.26) and (A.34), we get the L.H.S of Eq. (A.27) 

iT 2 iT 3 

L.H.S = m, + (5fO,,3B,^3 W, + ^O.y.zB^.^l^^ + ^<1),>,iB3.^ W',. (A.36) 
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By Eqs. (A.34), (A.23) and (A.3), we gain the R.H.S of Eq. (A.27) 



R.H.S = (1 - Si)mi + "VijWj + StAir,l^rj,3Wj + St%r,2^rj,2Wj + St^Air,3^rj,lWj + 0(St^). (A.37) 



So, we have 



1 Sf o ^t^ ^ O ^ A 

mt=-(-St^ir,3Brj,3-^^iraBlj,2-^^ir^Bl^^^ (A.38) 

Let 

and for J + 1 < / < A/^ 

1 Sf ^ SP ^ 9 ^ 

^ijA = -(-Smr,3Brj,3 " ^<^/r,2^,y,2 " ^^/r,l^'y,i + 'I'/j + StAir,l^rj,3 + ^rA,v,20,y,2 + ^^ A.v^O.yj), (A.39) 

and we restrict the truncated error of Eq. (A.39) equal to OiSt"^) and get 

3 

0,,;4 = ^^^^C0eff(0,,;4,^^/). 

So, we have 

mj = ^ijAWj + 0(St^). (A.40) 

Appendix B. The coefficient matrices of the higher-order L-NSE with the zero-mean flow 

(A) The coefficients of Sf are given by the following matrix 

~ Is ^ V -^ '^ y ) ~Ts y V ^ '^ ^y ) 

^ kx {-kx^ - ^/ + 3 ky^CTe^ + 3 k^^CTe^ + 3 ky^CTy^ + 3 /:;c V^^) 

-:^ky [-k/ -ky^ -\-3 ky^CTe^ + 3 y^X V,2 + 3 y^3; Vv^ + 3 y^;^ Vy^) 

(B.l) 

(B) The coefficients of Sp are given by the following matrix 

Coeff(5,5[l, 1], SP) = — ky^cTe + — k/ky^cTe + — /:xV, + — ky^cTy + — ^/o-v + — /://:3;Vv (B.2) 
Coeff(5,5[l,2],^^^) = 0,Coeff(5,5[l,3],^r^) = 0, Coeff (5,5 [2, 1],^^^) = (B.3) 

C0eff(5,5[2,2],^r^) = -| k/ky^CTyCTeCrq + ^ ^/cr^O-^O-y + ^ k^^ky^CTy^CTq - 1/9 kjc^ky^CTy^CTe - ^ k/ky^CTy- 

^k^k^rr -^k "^(T CT^ + ^k ^rr^rr - -L k "^rr - -L k "^rr - ^ k "^rr rr rr - ^ k ^k ^rr rr rr - 
108 ^x '^y ^-^ e (^ '^x ^ v ^ e ~ ^ '^x ^ e ^ q j^qS -^ ^ 108 -^^ 9'^x'^e'^e'^^ (^ '^x '^y ^ e ^ e ^ q 

36 y ^ 9 "^y ^ y ^ g'^x'^y^v^e^q g'^x^v^e^q^^'^x'^y^e^q g i^x i^y ^ e ^ v ^ i<^ i^x ^ v ^ q 

9 l^x ^v ^ e ~ g '^x % ^y "I" 9 % ^v ^q "*" YS -^ v ^ ^ 

(B.4) 

Coeff (5,5 [2, 3], ^^^) = \ kjc ky^O-y^O-q - ^k^ ky^O-y Cre(Tq-\ kjc ky^CTe (TeO-q+l k^^ky O-y^CTq + | k^^ky (Je (T q (Ty- 
9 i^x % (^e ^£ ^q "I" 6 ^x % <^e <^^ "I" 9 ^x % ^v ^e ^^ "I" g ^x % <^e ^q "I" 9 ^x % ^e ^v (^q ~ 54 ^x % ^e ~ 
Yog % /Cy (Ty — Yog % /Cy (Ty — 27 % ^3; CT^ ~ 9 % % <^^ <^y + 9 % % <^y + 9 % % <^y ~ 9 % % <^v <^e + 

To /C;^; /Cy (T £ g Kx Ky (Ty (T g g Kx Ky (Ty (T g 

(B.5) 

Coeff(5,5[3, 1],^^^) = (B.6) 
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j_ 

"27 ' 



CoefF(5,5[3, 2], SP) = j^ k^^ky Cr^ - | k^^ky Cr^ CrqCry+\ k^ ky^CTy Cre(Tq-\ k^ ky^CTe (TeCTq- I k^^ky (Te (Te (Tq + 
g % % O'y (Tq — g Kx Ky (Ty (Tg — 27 % % O'e ~ JqS -^ y ^^ ~ T08 -^ y ^^ ~ 54 -^ % ^^ "*" 9 '^x '^y ^v + 
9 ^x % ^e ^v ^^ "I" 9 ^x % ^^ ^e ^^ ~ g f^x % ^v ^e "I" g ^x % ^v ^^ ~ 9 ^x % ^e ^v "I" g ^x % ^e ^^ "I" 

9 A^x % ^v ~ 9 ^x % ^y ^e ~^ 6 ^x % ^e ^^ 

(B.7) 



108 -^ V ^ 



Coeff (fi 5[3, 3], <5f3) = - 1 k/ky^o-y o;o-q + ^ ^/^y V^V^ - i k/ky^o-y^cTe - ^ ^x^^yVv - ^ ^/fc 

35 ^x ^v ~ J^ / -^ ^x % <^e ^e ^q ~ JoS 3^ ^^ "*" 9 % ^^ ^^ ^^ "I" "9" ^x % <^y ^e ^^ "*" 6 3^ ^^ ^ ^ ~ 

1 )L 4 2 _ 1 i^ 4 2 _ 1 ^ 4 3 _ J_ ^ 4 1 r 2r 2 2 _ 1 -l 2r 2 2 . 2 r 4 2 

g Ky u g u y g K,y u y u Q g i^x ^ V ^QS "^J ^ e ^ (^ i^x i^j ^ 6 ^ q g'^x'^y^e^v^g'^x^v^q 

Q Kx K,y (Ty + Tq /Cy (Xy (X^ q /Cy (Ty (T g (Tg q Ky (T £ (T€ (T n + "Tq '^X '^J (T £ 

(C) The coefficient matrix of dt"^ are given by the following matrix 

Coeff(5,5[l,l],^^'^) = (B.9) 

Coeff (5 5 [1, 2], dt"^) = -^i'kx{2{) ^/^y V^ cr, + 10 kx^(Ty cr, - 90 ^y V, (Tq kx^ + 10 kJ^(Ty cr, + 4 ^/+ 

^ ^ (B.IO) 

7 y^y"^ - 30 y^y Vy (Tq + 60 y^;, Ve O-^ ky^ + 14 kx^ky^ - 30 kx^ky^CTy (Tq) 

Coeff(5,5[l, 3], dt^) = -^ i • ky (60 ^/o-^ (Tq ky^ + M^/^y^ - 90 ^yV, (Tq kx^ + 20 ^/^yV^ 0-, + 7 ^Z- 

(B . 1 1 ) 

30 k^ky^CFy (Tq - 30 ^jc^CTy CT^ + 10 ^j^^CTy (Tg + 10 /:y^(Ty (T^ + 4 ^y^ j 
Coeff(5,5[2, 1], ^^4) = i . (^4^/,^^4 ^ ^^ ^^3^2 ^ ^^ ^^5) ^g^^2) 

Coeff(5,5[2, 2], ^^'^) = 0, Coeff(5,5[2, 3], 6t^) = (B.13) 

Coeff(5,5[3, 1], 6t^) = i . {Uik% + &3/^'/^y + ^0,5/:') (B.14) 

/ _1^2 19 ^2 1_4 5__2_,7_3^,1^3^,1^2_2 2_.2_2 1^^3 1^^, 

7227 2725 1 17 7 1 5 2 

^(Ty (T^ + 2j(Te(Tq(Ty + -^jCTgCrq (Ty + jogCr^CTy + jogCTgCr^ — y^^CTyCrg — -^^(TyCTq + YgCT^CT^ — 2j(Te(Tq(Ty — 

52 12 1 27 212 14 13 31 21 3 

27^6 ^v(Tq ~ 27^^ ^q^e ~ 2j^e^q(Te "I" 2j^v(TQ(Te ~ ^j^q ^e(Te ~ 27^^ "*" 1620 ~ 648^^ ~ 27^^^^ 

(B.15) 

/ _13,1_2 4_2 2_4 4__2_,8_3^,2_3^,1^2_2 4_2_2 2__3 
^2,3 ~ sTO 54^ ~ 81^ ~ 27^ ~ fl ^ 'i ^ Tf^ ^ 9^ ^ 9^ ^ ~ 27^ ^ ~ fj^v^Te — 

13 4 2 2 11 2 11 2 1 1 25 13 1 

YQ^CTgCTq + 27 ^v <^^ "■" 27^^^^^^ 27^^^^ ^^ 27^^^^ 54^e<^e ~ 324 ^y^e ~ T08^^^^ Ts^^^^" 

4 242 22 2 2 11 222 24 19 22 3 

2j(T^(Tq(Ty — 27^6 ^y(Tq ~ 27^^ ^q^e ~ 27 ^^^^^^ "*" 27^^^^^^ ~ 27^^ ^ e(T e ~ 27^^ ~ yiA^^ ~ 27^^^^ 



fo,5 = -^ ^6^ ^y^ + ^ ^yV^2 + ^ Cr^ ^^^^v " ^ ^^'^^y^ + ^ ^e^^y ^^ + 1/9 0"^ V^ - ^ O"/ + ^ CTy V^ 
1/10^2^ 2 2___2 l^^^l^^ 2__2^^1^^^2 1^^^: 



27 CTy^ +1/18 (T^(Tq - i^(Te(Tq (Ty^ - j^ (T^ (Ty -\- j^ (Te (T^ - 2j ^e (Tq^CTy + 7^ (T^CTq (Ty^ - 7^ (T^CTq (Tg 

12 122 21 3112 525 

27 ^q ^e ^ e ~ 543 ^e ~ 27 ^^ ^^ ^^ ~ 27 ^^ ^^ 1620 ~ 27 ^^ ^^ ^^ 648 ^^ 162 ^^ ^^ 

(B.17) 
Coeff(5,5[3,2],^^'^) = 0,Coeff(5,5[3,3],^^'^) = (B.18) 

Appendix C. The coefficient matrices of tfie fiigfier-order L-NSE witfi tfie uniform fiow 

Coeff(5,2[l, 1], ^0 = 0, Coeff(5,2[l, 2], 6i) = 0, Coeff(5,2[l, 3], St) = (C.l) 



Coeff(5,2[2, list) = -Iky^UV^cTy + I ky^UcTy -3kxky U^VcTy + kx^UV^cTy + kxky V^cTy + I kx^UcTy- 
kx^U^cTy - kx^UV^cTe + I /:x^^^. - kx^U^cTe - kx ky U^VcTe + \kxky V(Te - kx ky y V, 
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(C.2) 



Coeff(5,2[2, 2], St) = -\k^^(Ty + ^/vV^ - \ ky^cr^ - \ k^W^cTy + | k^U^dy + ?>k^ky UVcry-\- 

k, ky UV(Te + I k,^V^(Te " | /:/cr, + I /:/^V, 
Coeff(5,2[2, 3], ^0 = 2 ky^UVcTy +\k^ky U^cTy -Ik^ky V^o-y - k/UVcTy + \k^ky U^cTe + k^^UVcTe- 

^k / 



^x % ^e "I" 2 -^ y ^^ 



COeff(5,2[3, 1], ^0 = -3 ky kj, UV^CTy - ky^V^CTy + I ^;c^ VO-v + ky^U^VCTy + ^3; ^;, U^ (Ty + ^ ^/VCTv- 

2 k/VU^cTy - ky^V^o-e - ky kj, UV^o-e + \ ky^Vo-e - ky kj, U^cr, + \kykx Ucr, - ky^U^Vcr, 
Coeff(5,2[3, 2], St) = -ky^UVcTy -jk^ky U^cTy +lkj,ky V^cTy + 2 k^^UVcTy -\kxkycre + lkj, ky U^cre+ 

Ik^kyV^O-e'tky^UVo-e 

Coeff(5,2[3, 3], St) = -\ k^^cTy -\-3kxky UVcTy - | ky^cTy + k^^U^cTy - \ ky^U^o-y + I ky^V^o-y- 

\ ky^o-e + kj, ky UVcr, + I y^/yV, + ^ y^/^V, 



(C.3) 
(C.4) 

(C.5) 

(C.6) 

(C.7) 



Appendix D. The optimized values of free parameters 

Considering n = 5 and w^^ = in Eqs. (37) and (59), cr^ = 0.0025 and CTy = 0.0025, the analytic expressions of 
the problems (64) and (67) are given by 

F%E) = -0.92763125500-^ + 0.1201533868 cr, + 7.203744088 cr/ + 129.3388220 o-,o-^ - 0.4322246453 0-^0-/+ 
0.01296673935 o-^V,^ - 0.00007699001495 cr/cr, + 11 3.4657225 o-,V/ + 43.21814228 0-^0-,^ + 
0.006483369681 cr^V/ + 0.01 836896999 cr^V, - 2. 155399630 cr^V, - 1.728866164 o-^V/+ 
0.006483369681 cr/o"/ - 0.5121252003 o-^o-,^ - 1.296718509 o-^V,2+ 
0.1080561613 o-,^ + 0.009393468450 cr^^ + 16.02085679 o-,^ - 0.00006398720339 cr^^ + 57.80191535+ 



0.00000023046353 16 o-/ 



(D.l) 



F'(E) = 0.007225929590 0-/0-^2 + 8.028810655 cr/ + 0.0000002866536304 0-^2 _ 0.1773049091 cr,- 

0.000035251935840-^ - 0.4817286393 o-^o"/ - 0.00008580791386 o-^V,+ (^ 2) 

0.008179411071 o-,o-^ + 0.001159746366 
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